

/*
csolve(pointer matrix f, numeric matrix x, real scalar crit, real scalar itmax)
f should be written so that any parametric arguments are packed in to x,
and so that if presented with a matrix x, it produces a return value of
same dimension of x.  The number of rows in x and f(x) are always the
same.  The number of columns is the number of different input arguments
at which f is to be evaluated.

crit:     if the sum of absolute values that f returns is less than this,
           the equation is solved.
itmax:    the solver stops when this number of iterations is reached, with rc=4
rc:       0 means normal solution, 1 and 3 mean no solution despite extremely fine adjustments
           in step length (very likely a numerical problem, or a discontinuity). 4 means itmax
           termination.
This is a translation by Matt Weinberg of matlab code originally written by Christopher Sims.  Matlab code available at
www.princeton.edu/~sims.  Any errors are the responsibility of Matt Weinberg.
mweinber@princeton.edu.

*/

cap mata: mata drop csolve()
mata
    numeric matrix csolve(pointer(function) scalar f, numeric matrix userx, real scalar crit, real scalar itmax,|a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15)
    {
	x=userx
	real scalar delta, alpha, verbose, done, itct, randomize, factor, shrink, subDone, lambda, lambdamin
	real matrix tvec
	numeric matrix f0, af0, af00, grad, dx0, dx, z, af, fmin, xmin
	transmorphic p
    p=callf_setup(f,args()-4,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15)    
    
    
    /*differencing interval for numerical gradient*/
    delta=1e-6
    /*************alfa*************************/
    /*tolerance on rate of descent*/
    alpha=1e-3
    /****************verbose****************/
    verbose=1 /*if this is set to zero, screen output is suppressed*/

    /*need to rewrite code for analytic derivates*/

    nv=max((rows(x),cols(x)));
    tvec=delta*I(nv);
    done=0;
    f0=callf(p,x)
    af0=colsum(abs(f0))
    af00=af0
    itct=0

    while(!done){ 
        if(itct>3 & af00-af0<crit*max((1,af0)) & mod(itct,2)==1){
            randomize=1
        }
        else {
            grad = (callf(p,x*J(1,nv,1)+tvec)-(f0)*(J(1,nv,1)))/delta
            if(isreal(grad)){
                if(1/cond(grad)<1e-12){
                    grad=grad+tvec
                }
                dx0=-luinv(grad)*f0
                randomize=0
            }
            else {
                if(verbose) display("gradient imaginary")
                randomize=1
            }
        }
        if(randomize){
            if(verbose) printf("\n Random Search")
            dx0=norm(x):/invnormal(uniform(rows(x),cols(x)))
        }
        lambda=1
        lambdamin=1
        fmin=f0
        xmin=x
        afmin=af0
        dxSize=norm(dx0)
        factor=.6
        shrink=1
        subDone=0
        while(!subDone){
            dx=lambda*dx0
            z=callf(p,x+dx)
            af=colsum(abs(z))
            if(af<afmin) {
                afmin=af
                fmin=z
                lambdamin=lambda
                xmin=x+dx
            }
            if(((lambda>0) & (af0-af<alpha*lambda*af0))| ((lambda<0) & (af0-af<0))){
                if(!shrink){
                    factor=factor^.6;
                    shrink=1;
                }
                if(abs(lambda*(1-factor))*dxSize>.1*delta){
                    lambda=factor*lambda
                }
                else if((lambda>0) & (factor==.6)) lambda=-.3
                else {
                    subDone=1
                    if(lambda>0) {
                        if(factor==.6) {
                            rc=2
                        }
                        else {
                            rc=1
                        }
                    }
			else {
				rc=3
			}
                }
            }
            else if((lambda>0) & (af-af0>(1-alpha)*lambda*af0)){
                if(shrink) {
                    factor=factor^.6
                    shrink=0
                }
                lambda=lambda/factor
            }
            else {
                subDone=1
                rc=0
            }
        }/*while(!subdone)*/
        itct=itct+1
        if(verbose) {
            printf("\nitct %f, af %g, lambda %g, rc %g\n",itct,afmin,lambdamin,rc) 
	    for (k=1; k<=rows(xmin); k++) {
            	printf("x %g\n", xmin[k])
	    }
            for (k=1; k<=rows(fmin); k++) {
		printf("fmin %g\n", fmin[k])
	    }
        }
        x=xmin
        f0=fmin
        af00=af0
        af0=afmin
        if(itct>=itmax) {
            done=1
            rc=4
        }
        else if(af0<crit) {
            done=1
            rc=0
        }


    }
    return(x\rc)
}

end
